/*!  \file  BiSVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/02/2013 04:23:54 PM
 *   \copyright GNU Public License.
 */
#include "BiSVD.h"
#include "TriEigen.h"
#include "Matrix.h"
#include <stdio.h>
#include "MatrixGenerator.h"
#include "commonSVD.h"
#include "ThreadPool.h"
BiSVD::BiSVD(const char * name_): Implementation(name_, BiSVDGroup), eigen(0),
		aUpper(false), aLower(false), aURowMajor(false), aUColMajor(false), 
		aVRowMajor(false), aVColMajor(false), VColVec(true){
		canBeParallel = false;
		setCompileTime(__DATE__, __TIME__);
		
};
void BiSVD::setTriEigen(TriEigen & implementation) { 
		eigen = &implementation; 
		setCompileTime(eigen->getCompileTime());
		canBeParallel = eigen->getCanBeParallel();
		nThreads      = eigen->getNThreads();
};

bool BiSVD::setNumberOfThreads(int nThreads_) {
	if (eigen) { 
		bool res = eigen->setNumberOfThreads(nThreads_);
		nThreads = res ? nThreads_ : nThreads;
		return res;
	}
	return false;
};
bool BiSVD::checkInput(const Matrix & A, Matrix & U,
 	 						  const Matrix & S, Matrix & V){

	assert(aUpper     || aLower);
	assert(aUColMajor || aURowMajor);
	assert(aVColMajor || aVRowMajor);
	
	bool res  =   A.rows() ==  U.rows() &&   U.rows() == U.cols() && 
	U.cols() == S.rows() && S.cols() == V.rows() &&  V.cols() == A.cols();
	
	if (!res) {
		cerr << "[BiSVD::checkInput] Error dimentions don't match" << endl;
		return res;
	}
	res &= A.getData() != NULL &&  S.getData() != NULL && 
		    U.getData() != NULL &&  V.getData() != NULL;

	if (!res) {
		cerr << "[BiSVD::checkInput] Error Matrices are not allocated" << endl;
		return res;
	}
	if ( U.isRowMajor() && !aURowMajor && U.isFullMatrix()) U.transpose();	
	if ( U.isColMajor() && !aUColMajor && U.isFullMatrix()) U.transpose();	
	if ( V.isRowMajor() && !aVRowMajor && V.isFullMatrix()) V.transpose();	
	if ( V.isColMajor() && !aVColMajor && V.isFullMatrix()) V.transpose();	
	if (U.isColMajor() && !aUColMajor) {
		cerr << "[BiSVD::checkInput] U should be in ColMajor format" << endl;
		return false;
	}
	if (U.isRowMajor() && !aURowMajor) {
		cerr << "[BiSVD::checkInput] U should be in RowMajor format" << endl;
		return false;
	}
	if (V.isColMajor() && !aVColMajor) {
		cerr << "[BiSVD::checkInput] V should be in ColMajor format" << endl;
		return false;
	}
	if (V.isRowMajor() && !aVRowMajor) {
		cerr << "[BiSVD::checkInput] V should be in RowMajor format" << endl;
		return false;
	}

	return true;
}
void BiSVD::saveCompute(const Matrix & B, Matrix & U, Matrix & S, Matrix & V) {
	reset();
	// check if input task is VALID
	success = checkInput(B, U, S, V);
	assert(B.isFullMatrix()); // for the moment
	if (!success) return;

	if ((!aUpper && B.isUpper()) || (!aLower && B.isLower())) {
		Matrix BT(B, Transpose);
		S.transpose();	
		saveCompute(BT, V, S, U);
		U.transpose();	
		S.transpose();	
		V.transpose();
		assert(checkSVDDimentions(B, U, S, V));
		return;
	} else { // B is OK :) 

		const double * RS d = B.getData();
		const double * RS e = B.getOffDiagonal();
		const double dLimit = 1e-13;
		const double eLimit = 1e-13;
		const int         l = B.minDim();

		int dLow[l]; int nDLow = 0;
		int eLow[l]; int nELow = 0;

		for (int i = 0; i < B.minDim() - 1; i++) {
			dLow[nDLow] = i;
			nDLow += fabs(d[i]) < dLimit;
			eLow[nELow] = i;
			nELow += fabs(e[i]) < eLimit;
		}
		dLow[nDLow] = l - 1;
		nDLow += fabs(d[l-1]) < dLimit;
		if ( nDLow > 0 || nELow > 0) { // we can split the problem
//			cout << "Spliting problem nDLow = " << nDLow << endl;
			if ( dLow[0] == 0 ) { // special case 
				if (nDLow + nELow == 1) {
					Matrix L   (l, l, Orthogonal);// B_new = L * B_org 
					L.makeIdentity();
					Matrix Bnew(B,  CopyWithData);// B_new = U * S * V;
					Bnew.deflateLeft(L);          // B_org = L.T() * U * S * V;

					assert((Bnew  - L * B).nrm2() < 1e-13);
					V.fill(0);
					Matrix BSub(Bnew); Matrix USub(U), VSub(V), SSub(S);
		
					USub.setView (1, 1, l-1, l-1);
					VSub.setView (1, 1, l-1, l-1);
 					SSub.setView (1, 1, l-1, l-1);
					BSub.setView (1, 1, l-1, l-1);

#ifndef NDEBUG	
					Matrix::Order orderU = U.getOrder(), orderV = V.getOrder();
#endif
					compute(BSub, USub, SSub, VSub);
					assert(orderU == U.getOrder());
					assert(orderV == V.getOrder());

//					double subNorm = VColVec ? 
//						( BSub - USub * SSub * VSub    ).nrm2() :
//						( BSub - USub * SSub * VSub.T()).nrm2();
//					assert(subNorm < 1e-13);

//					cout << "Sub  Norm = " << subNorm << endl;
//					cout << endl;
//					cout << "USub = " <<  USub << endl;
//					cout << endl << "VSub = " <<  VSub << endl;


					VSub.setView(0, 1, 1, l-1); VSub.fill(0);
					VSub.setView(1, 0, l-1, 1); VSub.fill(0);
					VSub.setView(0, 0, 1,   1); VSub.fill(1);

					USub.setView(0, 1, 1, l-1); USub.fill(0);
					USub.setView(1, 0, l-1, 1); USub.fill(0);
					USub.setView(0, 0, 1  , 1); USub.fill(1);
					SSub.setView(0, 0, 1  , 1); SSub.fill(0);
					U = L.T() * U;

//					double fullNorm = VColVec ? 
//						( B - U * S * V    ).nrm2() :
//						( B - U * S * V.T()).nrm2();
//					assert(fullNorm < 1e-13 * B.nrm2() * l);
				}
			}
	
		} else {
			if (B.rows() == B.cols())
				compute(B, U, S, V);
			else {
				assert(B.isUpper());
				U.fill(0); V.fill(0); S.fill(0);
				int m = B.minDim();
				int n = B.cols() > B.rows() ? B.cols() : B.rows();
				U.setView( 0, 0, m, m);
				V.setView( 0, 0, m, m);
				Matrix Bsub(B);
				Bsub.setView( 0, 0, m, m);
				S.setView( 0, 0, m, m);
				compute(Bsub, U, S, V);

				U.resetView();
				S.resetView();
				V.resetView();

				Matrix & big = U.rows() > V.cols() ? U : V;

				big.setView(0, m, m, n - m); big.fill(0);
				big.setView(m, 0, n - m, m); big.fill(0);
				big.setView(m, m, n - m, n - m); big.makeIdentity();
				big.resetView();

			}
		}
	}
	if (!VColVec) { 
		V.transpose();
	}
	assert(checkSVDDimentions(B, U, S, V));
}
inline double sqr(const double a) { return a * a; };

void BiSVD::computeNormal(const Matrix & B, Matrix & U, 
		Matrix & S, Matrix & VT ) {
	reset();
	unsigned int n = B.rows();
	const double * bd = B.getData();
	const double * be = B.getOffDiagonal();
	Matrix TU(n, n, TriDiagonal);
	Matrix TV(n, n, TriDiagonal);
	extraMemory = sizeof(double)*(TU.getStorageSize() + TV.getStorageSize());	
	double * du = TU.getModData();
	double * eu = TU.getModOffDiagonal(); 
	double * dv = TV.getModData();
	double * ev = TV.getModOffDiagonal();
	// compute B * B.T() 
	dv[0] = sqr(bd[0]);
	for (unsigned int i = 0; i < n-1; i++){
		du[i] = sqr(bd[i]) + sqr(be[i]);
		eu[i] = bd[i + 1] * be[i];
		dv[i+1] = sqr(bd[i+1]) + sqr(be[i]);
		ev[i] = bd[i] * be[i];
	}
	du[n-1] = sqr(bd[n-1]);
#ifndef NDEBUG
	double tol = 1e-12;
#endif
	eigen->compute(TU, S,  U);	assert((TU *  U -  U * S).nrm2() < tol);
	success = eigen->computaionSuccessed();
	eigen->compute(TV, S, VT); assert((TV * VT - VT * S).nrm2() < tol);
	assert(VT.isColMajor());
	extraMemory += eigen->getMaxUsedMemory();
	maxMemory = maxMemory < extraMemory ? extraMemory : maxMemory;
	success &= eigen->computaionSuccessed();
	double * s = S.getModData();
	for (unsigned int i = 0; i < n; i++) s[i] = sqrt(s[i]);
	// NOTE THAT if we multiply any column in U by -1, we preserve all 
	// properties of the eigen values decomposition
	//	Matrix TMP1 = U * S ; 
	//	Matrix TMP2 = B * VT;
	//	Matrix tmp = U.T() * B * VT;
	// tmp(i,i) = U'(i,j) * B(j,k) * VT(k, i)
	
	int nChanges = 0;
	int change[n];
	assert(B.isUpper());
	for (int i = 0; i < B.rows(); i++) {
		double res = 0;
		for (int j = 0; j < B.rows() - 1 ; j++) {
			res += U(j,i) * bd[j] * VT(j,i) + U(j,i) * be[j] * VT(j+1, i);
		}
		int j = B.rows() - 1;
		res += U(j,i) * bd[j] * VT(j,i);
		if ( res < 0 )	change[nChanges++] = i;
	}
	for (int i = 0; i < nChanges; i++) 	
		VT.changeColSign(change[i]);
	U.flipCols();	S.flipRows();	VT.flipCols(); 
	
//	assert(( B - U * S * VT.T()).nrm2() < B.nrm2() * B.minDim() * 1e-15 * 100);
}
void BiSVD::computeNormal2(const Matrix & B, Matrix & U, 
		Matrix & S, Matrix & V ) {
	reset();	
	unsigned int n = B.rows();
	const double * RS bd = B.getData();
	const double * RS be = bd + n;
	Matrix TU(n, n, TriDiagonal);
	extraMemory = sizeof(double) * TU.getStorageSize();	
	double * RS du = TU.getModData();
	double * RS eu = TU.getModOffDiagonal(); 
	for (unsigned int i = 0; i < n-1; i++){
		du[i] = sqr(bd[i]) + sqr(be[i]);
		eu[i] = bd[i + 1] * be[i];
	}
	du[n-1] = sqr(bd[n-1]);
	double tol = 1e-12;
	eigen->compute(TU, S,  U);	assert((TU *  U -  U * S).nrm2() < tol);
	success = eigen->computaionSuccessed();
	extraMemory += eigen->getMaxUsedMemory();
	maxMemory = maxMemory < extraMemory ? extraMemory : maxMemory;
	double * s = S.getModData();
	for (unsigned int i = 0; i < n; i++) s[i] = sqrt(s[i]);
	U.flipCols();
	S.flipRows();
	double * RS v = V.getModData();
	const double * RS u = U.getData();
	for (unsigned int i = 0; i < n; i++) { // compute V
		B.solve( u + U.offset(0,i), U.isColMajor() ? 1 : U.lda(),
				   v + V.offset(0,i), V.isColMajor() ? 1 : V.lda(),
					s[i]);
	}

/*	if (( B - U * S * V    ).nrm2() < B.nrm2() * B.minDim() * 1e-15)
		cout << "BRAVO 1" << endl;
	if (( B - U * S * V.T()).nrm2() < B.nrm2() * B.minDim() * 1e-15)
		cout << "BRAVO 2" << endl;
	if (VColVec) 
		assert(( B - U * S * V    ).nrm2() < B.nrm2() * B.minDim() * 1e-15);
	else
		assert(( B - U * S * V.T()).nrm2() < B.nrm2() * B.minDim() * 1e-15);	
*/
}
void BiSVD::computeGK(const Matrix & B, Matrix & U, 
		Matrix & S, Matrix & V ) {
	bool dimentionMatch = checkSVDDimentions(B, U, S, V);
	bool haveMem = checkSVDMemory(B, U, S, V);
	if (!eigen || !haveMem || !dimentionMatch) return;
	int n = B.rows();
	const double * RS bd = B.getData();
	const double * RS be = bd + n;
	Matrix GK(2 * n, 2 * n, TriDiagonal);
	double * RS d = GK.getModData();
	double * RS e = GK.getModOffDiagonal();
	for (int i = 0; i < 2 * n - 1; i++) {
		d[i] = 0.0;
		e[i] = i & 1 ? be[ i >> 1] : bd[ i >> 1];
	}
	d[2 * n - 1] = 0;
	int r = eigen->canComputeRange ? n : 2 * n;

	Matrix U1(2 * n, r, Orthogonal);
	U1.makeColMajor();
	Matrix S1(r    , r, Diagonal  );

	if ( eigen->canComputeRange ) {
		eigen->computeRange(GK, S1, U1, 1 , n);	
	} else 
		eigen->compute(GK, S1, U1);	
	assert((GK * U1 - U1 * S1).nrm2() < 1e-12);

	const double s2 = sqrt(2);
	double * RS  s =  S.getModData();
	double * RS  u =  U.getModData();
	double * RS  v =  V.getModData();
	double * RS u1 = U1.getModData();

	assert(U1.isColMajor());
	assert(V.isColMajor());
	assert(U.isColMajor());
	if ( eigen->canComputeRange ) {
		for (int i = 0; i < n; i++) {
			s[i] = -S1.getData()[i]; 
			double * RS tmpU1 = u1 + U1.offset(0,i);
			double * RS tmpU  = u  + U.offset(0,i);
			double * RS tmpV  = v  + V.offset(0,i);
			for (int j = 0; j < n; j++) {
				*tmpU = *(tmpU1 + 2 * j + 1) * s2;
				*tmpV = *(tmpU1 + 2 * j + 0) * (-s2);	
				tmpU++; tmpV++;
			}
		}
	} else {
		for (int i = 2 * n - 1; i >= n; i--) {
			const int c = 2 * n - 1 - i;
			s[c] = S1.getData()[i]; 
			double * RS tmpU1 = u1 + U1.offset(0,i);
			double * RS tmpU  = u  + U.offset(0,c);
			double * RS tmpV  = v  + V.offset(0,c);
   
			for (int j = 0; j < n; j++) {
				*tmpU = *(tmpU1 + 2 * j + 1) * s2;
				*tmpV = *(tmpU1 + 2 * j + 0) * s2;	
				tmpU++; tmpV++;
			}
		}
	}
//	cout << B - U * S * V.T() << endl;
}
BiSVDBenchmark::BiSVDBenchmark(BiSVD & bisvd_):Benchmark(bisvd_) {
	bisvd = &bisvd_;
	set(baseName, "Benchmark BiSVD");
	const char * s1 = baseName;
	const char * s2 = bisvd->getName();
	char s3[strlen(s1) + strlen(s2) + 2];
	sprintf(s3, "%s %s", s1, s2);
	setName 				( s3);
	setDescription  	("Measures performance of registered bi svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	nRuns = 10;
	nConfigurations = 15;
	setConfigurationEachRun = true;
	generator = new MatrixGenerator();
	generator->setZeroProb  (0.0);
	Bi = U = V = S = NULL;
	setCompileTime(__DATE__, __TIME__);	
	setCompileTime(bisvd->getCompileTime());; 
// METRICS 
	addTimerMetrics();
	addMemoryMetrics();
	addSuccessMetrics();
	for (int i = 0; i < nSvdNorms; i++) svdNormMetrics[i] = 0;
	addSVDNormsMetrics();
}
BiSVDBenchmark::~BiSVDBenchmark() {
	if (Bi) delete Bi; Bi = 0;
   if (U ) delete U ; U  = 0;
   if (S ) delete S ; S  = 0;
   if (V ) delete V ; V  = 0;
	if (generator) delete generator; generator = 0;
//	for (int i = 0; i < nSvdNorms; i++) svdNormMetrics[i] = 0;

}
int BiSVDBenchmark::getProblemSize(const int configuration) const {
	const int problemSize[] = { 10, 25, 50, 75, 100, 150, 200, 300, 400, 500, 600, 700, 1000, 1500, 2000, 2500, 3000, 4000, 5000, 6000, 7000, 8000 };
	return problemSize[configuration];

};
void BiSVDBenchmark::setRunConfiguration(const int n) {
//#ifdef NDEBUG
	int N = getProblemSize(n);
	if (Bi) delete Bi; Bi = 0;	Bi = new Matrix(N, N, BiDiagonal);
	if (U ) delete U ; U  = 0;	U  = new Matrix(N, N, Orthogonal);
	if (U ) U->makeColMajor();
	if (S ) delete S ; S  = 0;	S  = new Matrix(N, N, Diagonal  );
	if (V ) delete V ; V  = 0;	V  = new Matrix(N, N, Orthogonal);
	
	if (Bi && U && S && V) {
		generator->generate(*Bi);
	}

}
void BiSVDBenchmark::prepareRun() {}
void BiSVDBenchmark::makeRun() {
	if (Bi && U && S && V && bisvd) { 
		bisvd->setNumberOfThreads(1);
		bisvd->saveCompute(*Bi, *U, *S, *V);
	}
}

void BiSVDBenchmark::addSVDNormsMetrics() {
	for (int i = 0; i < nSvdNorms; i++) 
	if (svdNormMetrics[i] == 0) {
		svdNormMetrics[i] = new Metric<double> (nRuns);	
		if (svdNormMetrics[i]) 
			addMetric(*(svdNormMetrics[i]), normNames[i]);
	}
} 

void BiSVDBenchmark::measureSVDNorms() {
	if (Bi && U && S && V ) {
		for (int i = 0; i < nSvdNorms; i++)
			if(!svdNormMetrics[i]) return;
		double svdNorms[nSvdNorms];
		getResidual(svdNorms[0], *Bi, *U, *S, *V);
		double tmp;
		getOrth    (svdNorms[1], *U); 	getOrth    (tmp, *V);
		svdNorms[1] = fmax(tmp, svdNorms[1]);
		for (int i = 0; i < nSvdNorms; i++)
			svdNormMetrics[i]->addMeasurement(svdNorms[i]);
	}

}
void BiSVDBenchmark::finalizeRun() {
	measureSVDNorms();
}

BiSVDScaling::BiSVDScaling(BiSVD & bisvd_) : Benchmark(bisvd_){
	bisvd = &bisvd_;
	set(baseName, "Scaling BiSVD");
	const char * s1 = baseName;
	const char * s2 = bisvd->getName();
	char s3[strlen(s1) + strlen(s2) + 2];
	sprintf(s3, "%s %s", s1, s2);
	setName 				( s3);
	setDescription  	("Measures scaling of registered bi svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	nRuns = 10;
	nConfigurations = ThreadPool::getOnlineCores();
	setConfigurationEachRun = true;
	generator = new MatrixGenerator();
	generator->setZeroProb  (0.0);
	Bi = U = V = S = NULL;
	int N = 500;
	Bi = new Matrix(N, N, BiDiagonal);
	B2 = new Matrix(N, N, BiDiagonal);
	U  = new Matrix(N, N, Orthogonal);
	S  = new Matrix(N, N, Diagonal  );
	V  = new Matrix(N, N, Orthogonal);
	if (U ) U->makeColMajor();
	if (Bi) generator->generate(*Bi);
	if (B2 && Bi) B2->copyWithData(*Bi);
	setCompileTime(__DATE__, __TIME__);	
	setCompileTime(bisvd->getCompileTime());; 
// METRICS 
	addTimerMetrics();
	addMemoryMetrics();
	addSuccessMetrics();
	for (int i = 0; i < nSvdNorms; i++) 
		svdNormMetrics[i] = 0;
	addSVDNormsMetrics();
	xlabel = "Number of cores";
}
BiSVDScaling::~BiSVDScaling() {
	if (Bi) delete Bi; Bi = 0;
	if (B2) delete B2; B2 = 0;
   if (U ) delete U ; U  = 0;
   if (S ) delete S ; S  = 0;
   if (V ) delete V ; V  = 0;
	if (generator) delete generator; generator = 0;
}
int  BiSVDScaling::getProblemSize(const int configuration) const {
	return configuration + 1;
}
void  BiSVDScaling::setRunConfiguration(const int n) {
	int threads = getProblemSize(n);
	if (bisvd && bisvd->getCanBeParallel()) {
		bisvd->setNumberOfThreads(threads);
	};
}
void  BiSVDScaling::prepareRun() {
	if (B2 && Bi) Bi->copyWithData(*B2);
}
void  BiSVDScaling::makeRun() {
	if (Bi && U && S && V && bisvd) { 
		bisvd->saveCompute(*Bi, *U, *S, *V);
	}
}
void  BiSVDScaling::finalizeRun() {
	measureSVDNorms();
}
void  BiSVDScaling::addSVDNormsMetrics() {
	for (int i = 0; i < nSvdNorms; i++) 
	if (svdNormMetrics[i] == 0) {
		svdNormMetrics[i] = new Metric<double> (nRuns);	
		if (svdNormMetrics[i]) 
			addMetric(*(svdNormMetrics[i]), normNames[i]);
	}
}
void  BiSVDScaling::measureSVDNorms() {
	if (B2 && U && S && V ) {
		for (int i = 0; i < nSvdNorms; i++)
			if(!svdNormMetrics[i]) return;
		double svdNorms[nSvdNorms];
		getResidual(svdNorms[0], *B2, *U, *S, *V);
		double tmp;
		getOrth    (svdNorms[1], *U); 	getOrth    (tmp, *V);
		svdNorms[1] = fmax(tmp, svdNorms[1]);
		for (int i = 0; i < nSvdNorms; i++)
			svdNormMetrics[i]->addMeasurement(svdNorms[i]);
	}
}

BiSVDUnitTest::BiSVDUnitTest(BiSVD & bisvd_): bisvd(bisvd_) {
	int l = strlen(bisvd.getName());
	char buf[l + 20]; 
	sprintf(buf, "Test BiSVD %s", bisvd.getName());
	setName 				(buf);
	setDescription  	("Test All registered bi svd implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
	setCompileTime(bisvd.getCompileTime());
};
bool BiSVDUnitTest::test1(const int m, const int n, const double * d, const double limit, 	const double * s_ref) const {
	bool res = true;
	const int  nNorms = 6;	double norm[nNorms];	
	// test with upper
	Matrix B(m, n, BiDiagonal, d); assert(B.isUpper());
   Matrix S(m, n, Diagonal,d);
	Matrix U(m, m, Orthogonal);
	Matrix V(n, n, Orthogonal);
	bisvd.saveCompute(B, U, S, V);
	res &= bisvd.computaionSuccessed();
	getNorms(norm, B, U, S, V );
	for(int i = 0; i < nNorms; i++) {
		res &= norm[i] < limit;
	}
	assert(res);
	if ( s_ref != NULL) {
	   Matrix Sigma_ref(m, n, Diagonal, s_ref);
		res &= (Sigma_ref - S).nrm2() < limit;
	}
	assert(res);
	// test with lower
	B.transpose();
	S.transpose();
	bisvd.saveCompute(B, V, S, U);
	
	res &= bisvd.computaionSuccessed(); 
	assert(res);
	getNorms(norm, B, V, S, U );
	for(int i = 0; i < nNorms; i++) {
		res &= norm[i] < limit;
	}
	if ( s_ref != NULL) {
	   Matrix Sigma_ref(n, m, Diagonal, s_ref);
		res &= (Sigma_ref - S).nrm2() < limit;
	}


	return res;
}
bool BiSVDUnitTest::test() const {
	bool res = true;	double limit = 1e-12;
	double d [] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
		1.0/8, 2.0/8, 3.0/8,	4.0/8, 5.0/8, 6.0/8, 7.0/8, 8.0/8, 9.0/8};
	double s_ref []       = {1.02891299546928341e+01, 
	9.00724709978768523e+00, 7.97107425087430066e+00, 6.97271747345694592e+00,
	5.97661256325545853e+00, 4.98057100904744665e+00, 3.98454752282448954e+00,
   2.98856032475197919e+00, 1.99267289366755218e+00, 9.97383294572933621e-01};
	res &= test1(10, 10, d, limit, s_ref);

	double d2[] =  {1, 2, -3, 1, -5, -6, 7, 8, 9, 10, 
		1.0/8, 2.0/8, 3.0/8,	4.0/8, -5.0/8, 6.0/8, -7.0/8, 8.0/8, 9.0/8};
	res &= test1(10, 10, d2, limit);
	res &= test1(6 , 5 , d2, limit * 1000);

	const double k = 1e-2;
	double dHard[] =  {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
		k * 1, k * 2, k * 3,	k * 4, k * 5, k * 6,	k * 7, k * 8, k * 9};
	
	res &= test1(10, 10, dHard, limit);
	double blow[] = {	-1.9972713886700525e+01, 9.7842422568977785e+00,  
							 1.1412294768584568e+01, 1.8900882963185474e+01,  
							 1.7502780876828847e+01,-5.3450974367118418e+00,
						   -5.2754827086988696e+00,-3.2518924437190715e+00, 
							-4.9932281621924197e+00  };
	res &= test1(5, 6, blow, limit * 10);

	return res;
	
	



/* 
	const unsigned int n2 = 26;
	MatrixGenerator generator;
	generator.setZeroProb  (0.0);
	Matrix B2 (n2, n2, BiDiagonal );
	generator.generate(B2);
   Matrix S2(n2, n2, Diagonal);
	Matrix U2(n2, n2, Orthogonal);
	Matrix V2(n2, n2, Orthogonal);
	bisvd.saveCompute(B2, U2, S2, V2);
	limit = B2.nrm2() * n2 * 1e-16 * 20;
	getNorms(norm, B2, U2, S2, V2 );
	for(int i = 0; i < nNorms; i++)
		res &= norm[i] < limit;
*/
/* double d3[] =  {0, 2, -3, 1, -5, -6, 7, 8, 9, 10, 
		1.0/8, 2.0/8, 3.0/8,	4.0/8, -5.0/8, 6.0/8, -7.0/8, 8.0/8, 9.0/8};
	B.initialize(d3);
	bisvd.saveCompute(B, U, S, V);
	res &= bisvd.computaionSuccessed();
	getNorms(norm, B, U, S, V );
	for(int i = 0; i < nNorms; i++) {
		res &= norm[i] < limit;
	}
*/
/*  
	double d4[] =  {0.1, 2, -3, 1, -5, -6, 7, 8, 9, 0, 
		1.0/8, 2.0/8, 3.0/8,	4.0/8, -5.0/8, 6.0/8, -7.0/8, 8.0/8, 9.0/8};
	B.initialize(d4);
	bisvd.saveCompute(B, U, S, V);
	res &= bisvd.computaionSuccessed();
	getNorms(norm, B, U, S, V );
	for(int i = 0; i < nNorms; i++)
		res &= norm[i] < limit;

	double d5[] =  {1, 2, -3, 1, 0, -6, 7, 8, 9, 10, 
		1.0/8, 2.0/8, 3.0/8,	4.0/8, -5.0/8, 6.0/8, -7.0/8, 8.0/8, 9.0/8};
	B.initialize(d5);
	bisvd.saveCompute(B, U, S, V);
	res &= bisvd.computaionSuccessed();
	getNorms(norm, B, U, S, V );
	for(int i = 0; i < nNorms; i++)
		res &= norm[i] < limit;
*/
	return res;
}
